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Abstract 

We introduce a new concept for generating optimal quadrature rules for splines. 
Given a target spline space where we aim to generate an optimal quadrature rule, 
we build an associated source space with known optimal quadrature and trans¬ 
fer the rule from the source space to the target one, preserving the number of 
quadrature points and therefore optimality. The quadrature nodes and weights 
are, considered as a higher-dimensional point, a zero of a particular system of 
polynomial equations. As the space is continuously deformed by modifying the 
source knot vector, the quadrature rule gets updated using polynomial homo¬ 
topy continuation. For example, starting with cubic splines with uniform 
knot sequences, we demonstrate the methodology by deriving the optimal rules 
for uniform cubic spline spaces where the rule was only conjectured hereto¬ 
fore. We validate our algorithm by showing that the resulting quadrature rule is 
independent of the path chosen between the target and the source knot vectors 
as well as the source rule chosen. 

Keywords: Gaussian quadrature, B-splines, well-constrained polynomial 
system, polynomial homotopy continuation 


1. Introduction 

Numerical integration of univariate functions is a fundamental mathemat¬ 
ical task which is a subroutine of many complex algorithms and is typically 
frequently invoked. Naturally, such an integration (or quadrature) rule must 
be as efficient as possible. We derive a new class of quadrature rules that are 
optimal in the sense that they require the minimal number of function’s evalu¬ 
ations. 
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A quadrature rule, or shortly a quadrature^ is an m-point rule, if m evalu¬ 
ations of a function / are needed to approximate its weighted integral over a 
closed interval [a, b] 

pb 'm 

/ w{x)f{x)dx = '^LOif{Ti) + Rmif), (1) 

i=l 

where re is a fixed non-negative weight funetion defined over [a,b]. The rule is 
required to be exaet, that is, Rm{f) = 0 for each element of a predefined linear 
function space C. The rule is said to be optimal if m is the minimal number of 
weights uji and nodes Ti, points at which / has to be evaluated. 

For the space of polynomials, the optimal rule is known to be the classical 
Gaussian quadrature [7] with the order of exactness 2m — 1, that is, only m 
evaluations are needed to exactly integrate any polynomial of degree at most 
2m — 1. Consider a sequence of polynomials {pQ,pi,... ,Pm^ • • •) that form an 
orthogonal basis with respect to the scalar product 

< f ,9 >= [ f{x)g{x)w{x)dx. (2) 

J a 

The quadrature points are the roots of the m-th orthogonal polynomial Pm 
which in the case when w{x) = 1 is the degree-m Legendre polynomial [13]. 

In this paper, we focus on piece-wise cubic polynomials, cubic splines, but 
the methodology is general and could be used for higher degrees as well. A 
univariate space of cubic splines is uniquely determined by its knot veetor. This 
knot vector is a non-decreasing sequence of real numbers called knots and the 
multiplicity of each knot determines the smoothness between the cubic pieces. 
To simplify the argument, we first study uniform knot vectors with all interior 
knots with uniform multiplieity. In particular, we investigate knot vectors with 
single interior knots which yield cubic splines and knot vectors with double 
interior knots which give cubic splines. For a more detailed introduction on 
splines, we refer the reader, e.g., to nn. 

The quadrature rules for splines have been studied since late 50 ’s [SllolIIl]. 
Firstly conjectured by Schoenberg [12], later proved by Micchelli and Pinkus |2|, 
the conditions of the existence and uniqueness of the optimal (Gaussian) rule 
have been derived. For spline spaces with maximum continuity (e.g., cubic 
splines), Micchelli and Pinkus |9| proved that there always exists a unique Gaus¬ 
sian quadrature rule. Their result, however, also reveals a nice phenomenon: 
the number of optimal nodes stays fixed as long as the number of interior knots 
stays constant. Therefore if one desires to derive a class of quadrature rules with 
the same number of nodes, spline spaces with higher continuity must have ad¬ 
equately more sub-intervals (elements) than spaces of lower continuity. This is 
natural since splines of lower continuity are limits of the higher continuous ones, 
when merging continuously two (or several) knots together, and their result is 
in agreement with this fact. 

For spaces with lower continuity, or when boundary constraints are involved, 
the rule is not guaranteed to be unique. Only spaces with fixed continuities 
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Figure 1: Spline space initialization flowchart 


were studied in [9]. To the best of our knowledge, the results on existence and 
uniqueness are not known for spaces with mixed continuities, we refer to these 
as knot vectors with mixed multiplicities. For example, for cubic splines, this 
includes to knot vectors with both single and double knots. In this work, we 
derive quadratures for this kind of mixed continuity spaces and show numerically 
that optimal quadratures exist. 

We use the observation that a spline space with multiple knots is a limit 
case of a spline space with the same number of single knots (when counting the 
multiplicities) when two, or several, knots merge. Abstracting this merging as 
a continuous transition between the two knot vectors of the same cardinality, 
the source and the target ones, the corresponding spline spaces continuously 
evolve from one into the other. The quadrature rule also depends continuously 
on the spline space since the quadrature rule can be seen as a zero (root) of 
a certain polynomial system and an infinitesimal change of the system does 
not significantly change the root. Based on these facts, we propose a new 
methodology that for a given (target) spline space S generates an associated 
source space S where the Gaussian quadrature rule is known. These spaces are 
defmed above knot vectors T jmd A, respectively, and the quadrature rule Q 
of S is numerically traced as X evolves into A, see Fig. 

The rest of the paper is organized as follows. Section briefly overviews 
the homotopy continuation of polynomial systems and Section summarizes a 
few basic properties of cubic splines. In Section we introduce a homotopy- 
continuation-based algorithm and discuss the results and the validity of the new 
quadratures obtained in Section Finally, Section summarizes our conclu¬ 
sions and describes future research directions. 


2. Homotopy continuation for polynomial systems 

Polynomial homotopy continuation (PHC) is a numerical scheme that solves 
polynomial systems of equations. This approach was introduced to solve the 
problem of movability of kinematic mechanisms [15] where the variables are the 
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free parameters of a certain mechanism tied together by a set of polynomial 
constraints. However, the unknowns and the constraints may relate to an arbi¬ 
trary problem. As our Gaussian quadrature rules are derived using a variation 
of this method, for the sake of completeness, we now briefly review the ideas 
used in homotopy continuation. For a detailed explanation, we refer the reader 
to the book m- 

Consider a well-constrained 2m x 2m polynomial system 

J'(x) = 0, X G C ( 3 ) 

where the domain is a hypercube in i.e., Q, = [x^xx] x • • • x [a;2m,3:2m]- 
Let K be an upper bound of the number of real roots of in and let us 
consider a simpler system *F(x) = 0 defined over some O with exactly K real 
roots {ri,..., yk}^ that is. 


J'(ri)=0, i = (4) 

For now, let us assume O = The roots {ri,..., ri^} are known, in fact they 
are chosen as an input since the system T is as simple as possible, see El- 
Consider now a continuous deformation of the known system into the desired 
one _ 

^(x) ^ ^(x). (5) 

Let us denote by H a one parameter family of systems created, as T is trans¬ 
formed into T . The parameter characterizing the transformation can be thought 
of as time or pseudo-time t. Thus 

?^(x, t)=0, X G tG[0,l], (6) 

where 'H(x, 0) is the artificially constructed system JF(x) = 0 , for which the 
roots are known, and ?^(x, 1) is the given system (§ which we seek to solve. 
Consider now the roots 0 as a function of time, ri(t),..., rx(t), and think of 
them as trajectories (curves) in As t runs over [0,1], some of these roots 

may vanish, which corresponds to the fact that 0 may have less real roots than 
LC, but no new trajectory may rise. The transformation is guaranteed not to 
introduce new roots. 

Typically, the source and the target systems are blended in a linear fashion, 
corresponding to the shortest path when deforming one polynomial system into 
another. However, one can consider different paths between the systems m- 


3. Transition between and cubic splines 


We recall several properties of spline basis functions. We consider a uniform 
knot vector 

JYn — (n — Xq , , ..., 1 5 1 5 (^) 
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Figure 2: Continuous transformation of cubic spline spaces. The uniform knot sequence with 
knots of multiplicity two is transformed to a uniform knot sequence consisting of single knots. 
Four corresponding basis functions of the source space an intermediate space, and the 

target space S ^2 shown. 


on the interval [a, 6], that is, each of the n — 1 interior knots has multiplicity two. 
We define h := =Xk— Xk-i for all /c = 1,..., n. Let us denote N := 2n — l 
and consider a uniform knot vector 


= {a = xq.xi, ... ,xn-i,xn = b), (8) 

each knot of multiplicity exactly one and define H := = Xk — Xk-i for all 

k = 1,..., TV. We denote by tts the space of polynomials of degree at most 3 
and define S'J the linear space of cubic splines over a uniform knot sequence 
as 

Sh = {/ e C^[a,b] : ^ 7r3,fc = (9) 

Similarly we denote by the linear space of cubic splines over a uniform knot 
sequence 

^3^2 = {/eC'2[a,6] e TTS,= (10) 


The dimension of both spaces is 2n + 2 = AT + 3. That is, the total number of 
interior knots is the same for both spaces, while the number of non-zero knot 
spans is different. 


Remark 1. The uniform knot sequences and are, in this work, taken as 
the source and the target knot sequences, respec^vely. Additionally, we consider 
the intermediate knot sequences generated as continuously evolves to T/v, 
see Fig. Given a particular knot vector pattern, which may contain both 
single and double knots, we also consider all the spaces of mixed continuities, 
not just ^ and (10). 


Consider now the space. Similarly to misiiii], we work with the non- 
normalized B-spline basis. To define the basis, we consider xq and Xn as double 
knots and extend our knot sequence with two extra double knots outside 
the interval [a, b] that we set to be 


X-i = xo — h and + h. 


( 11 ) 
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Figure 3: A target uniform knot vector with single knots for N = 5 elements and an 
associated source knot vector An with all knots of multiplicity two are shown. Both knot 
vector are of the same cardinality and have the same number of interior knots i = 4. 

see Fig. The choice of x_i and yields particular integral values for the 
first and last functions as discussed in ( |l3| . Denote by D = the basis 

of i where 

-^2/c — 1(0 ~ [^/c — 2? ^/c—2 5 1 5 ^/c — 1 5 ^/c] (• 

F^2/c(^) ~ — + 5 

where [.]/ stands for the divided difference and = niax(?i, 0 ) is the truncated 
power function. Direct computation gives 

/[5fc] = ifor fc = 3,4,...,2n, (12) 

where /[/] stands for the integral of / over the interval [a, b]. We work with non- 
normalized basis functions, and therefore the integrals above are independent 
on the knot sequence. With the choice made in 0 . direct integration gives 

I[Di] = I[D 2 n+ 2 ] = ^ and /[A] =/[A^+i] = A ( 13 ) 

Similarly to 0 , we extend the knot sequence of TV by two triplets of single 
knots as 


X-k=xo — kH and XN-\-k = /c = l,2,3. (14) 

and define D = the basis of S ^2 where 

^ki^^ ~ [^/c—4? ^k—St ^k — 2’! ^k — l'i ^/c] (• ^)+ = ‘2tI T 2. 

We obtain _ 

I[Dk]=I[Dk], A: = 4,...,2n-1 (15) 

and the six boundary integrals (three only due to symmetry) are computed 
directly by integrating Di, D 2 , and D 3 on [a, 6 ]. These integrals change during 
continuation and therefore have to be recomputed for various t. 

Now we consider a continuous transition between S^i and Fig. [2 

Since the transition of the spline spaces is governed by the transformation m 
the corresponding knot vectors, consider the following mapping 

Xn TV, (16) 
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Figure 4: Continuous transformation of knots sequences. The source knot sequence with 
double knots is changed over time and is transformed into the target uniform knot sequence 
with single knots, Ttv- The transformation is metaphorized as a path in 2u + 4 dimensional 
space of free knots. Two particular paths are shown: whereas one transformation changes 
always only a single knot (top), the geodesic transformation (bottom) generates the shortest 
path in the space of knots and consequently the shortest transformation from Xn into • 


including the six outer knots defined in (11) and 0, see Fig. Due to 
the fact that we work with non-normalized basis functions, Eqns. ( |15[ ) remain 
unchanged. The total number of knots is 2n + 6, but since the two boundary 
knots are constrained to stay fixed, there is 2n+4 free knots. The transformation 
can be conceptualized as a curve between and T/v, two points in see 

Fig-S There exist infinitely many paths connecting the source and target knot 
vectors. In particular, we analyze two specific knot transformations: the first one 
sequentially changes only one knot whilst all the others remain unchanged; the 
second one simultaneously spreads all free knots from the source position to the 
target one in a linear fashion. Such a transition can be seen as a diagonal straight 
line connecting and Tat, i.e., a geodesic path when under the Euclidean 
metric on the vector of free knots, while the first one corresponds to a path 
along the edges of a (2n + 4)-dimensional hypercube. 


Remark 2. Alternatively, one can extend and Tat by adding different six 
knots than those in 0 and ( p!4| ). Eor example, one can extend them in the 
open knot fashion by setting the multiplicity four to both boundary knots. In 
that case, the six (three due to symmetry) initial boundary integrals have to be 
computed accordingly. 
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4. Gaussian quadrature rules via homotopy continuation 

In this section, we use polynomial homotopy continuation (PHC) to generate 
Gaussian quadrature rules for families of cubic spline spaces above arbitrary 
knot sequences. 

The key observation is the following. _Assume you have an exact ar^ op¬ 
timal quadrature rule for a spline space S defined above a knot vector X and 
consider a continuous transformation of A as a function of time, i.e., X(t). As a 
quadrature rule is a solution of a particular polynomial system defined above A, 
the quadrature nodes and weights also change continuously. We use the PHC 
framework, to take an exact and optimal rule as the initial root of a certain 
polynomial system and trace this high-dimensional point (“follow the path of 
the root” in the PHC terminology, see El) while changing continuously X. 
This idea is in accordance with the result of Micchelli and Pinkus [9] which 
states that there exists an optimal quadrature rule with m nodes such that 


d -j- i “h 1 — ‘I'm 


(17) 


where d is the spline degree and i is the number of interior knots, counted 
including their multiplicities. Therefore m stays fixed as long as the number of 
interior knots remains constant. Homotopy continuation transfers the optimal 
quadrature rule from one space to another because the number of interior knots 
remains unchanged. 

In this paper, we want to derive optimal quadrature rules for spaces of 
cubic splines, so for uniform knot vectors we have S ^2 target spaces, see 

(10). The multiplicity of the interior knots is one and therefore Eq. 0 requires 


N to be odd. Further, we have n = and set S^i (|9|) as our source spaces 

For these spaces, unique quadrature rules that are explicit were derived in m- 
The framework we present is general and one can derive new rules by applying 
PHC to any other appropriate source space where an optimal quadrature rule 
is known, e.g., quintic splines with uniform knot sequences [3]. 


4-1. Gaussian quadrature formula 

We set our source space as see ^ and know, according to (|17|), that 
m = n + 1. The optimal source quadrature rule is of the form 


Q'alf] = 


nh Ti-\-\ 

/ /(i)di = i 

i=i 




(18) 


and the nodes and weights can be computed by ajecursion derived by Nikolov 
m, see Fig. Due to the equal dimensions of S^i and target rule 

requires the same number of nodes and therefore 




source space S^i 




target space S ^2 


•[r2,CJ2] • 


source rule Q 


[r2,(j02] 


target rule Q 


a Ti 


Xl T2 


X2 


Xl T2 


X4 b = X5 


Figure 5: For a desired (target) space S ^2 with N = 5 elements on [a, b], an associated source 

space S '3 3 ^, n = 3, is built (left). The source knot sequence (blue) is uniform with two interior 
knots Xl and X 2 i each of multiplicity two. The target kno t seq uence has four single knots, so 
the total number of interior knots is unchanged; z = 4 in Therefore both optimal rules 

require m = n + 1 nodes (green) and weights (red). The target rule (right) is derived via 
homotopy continuation, see Algorithm^ 


During the continuation, we transform the spline space S^i to 8^2 
accordingly the optimal rule Q ^ Q. Therefore Q, represented by the its nodes 
and weights, is a function oft. To simplify notation, if no ambiguity is imminent, 
we omit the time parameter and write instead of Ti{t). The source rule is 
Q = 2(0) and the target rule is Q = 2(1). We denote by 2 the rule (linear 
operator) and later in Section 4.3 use the symbol r for a 2m-dimensional point, 
but they both refer to the same set of optimal nodes and weights. Before we 
proceed to our homotopy continuation setting, we need to establish notation 
that unifies the source and the target quadrature rules. 

Let Ttv = {a = xo,xi,... ,Xiv-i,^Ar = 6) be a knot vector consisting of N 
subintervals, some of the subintervals may degenerate to zero length in the case 
of a double knot and let {n,... ,r^} be the quadrature points. We define a 
nodal pattern p of the quadrature rule 2 on T/v as an Wdimensional vector 
where the i-th coordinate specifies the number of quadrature nodes inside the 
i-th sub-interval. We say p is regular if no node coincides with a knot. 


Example 4.1. Consider the Gaussian quadrature rule of Nikolov m for § 
in the case when n is odd, n > 1. The rule says every sub-interval contains 
exactly one node except the middle one which contains two, see Fig[^ for n = 
3. Considering the multiplicities, the contains n sub-intervals of non-zero 
length and {n — 1) degenerated subintervals (double knots). The nodal pattern 
is then 


P = 


(1,0,...,1,0, 

n — 1 



0 , 1 ,..., 0 , 1 ) 

-V-^ 

n — 1 


( 20 ) 


and the number of nodes is exactly m = n 1. In the case when n is even, the 
middle node coincides with the middle knot. In such a case, we say the nodal 


pattern is irregular. We will discuss these patterns later in Section 4.3 


The nodal patterns describe the layout of the quadrature nodes with respect 
to the knots and, as long as the pattern does not change, the pattern of the 
polynomial system to solve remains unchanged. The crucial part of any optimal 
quadrature rule is to derive a correct nodal pattern. The result of Micchelli and 
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Pinkus [9] states the number of optimal nodes and also the knot span in which 
each particular node lies. However, one cannot conclude the nodal pattern for 
S ^2 from their result. In our approach once we know the nodal pattern of the 
source rule, then we know the initial polynomial system and its root. 


4 . 2 . Homotopic setting 

We are now ready to apply the PHC framework to the quadrature problem. 
The vector of unknowns consists of the quadrature nodes and weights 

X = (ti, . . . . . . ,Wm), X G 

our source polynomial system T expresses that the source rule Q exactly inte¬ 
grates the source basis D, that is, 

Ql{Di) = I[Di], i = l,...,2m (21) 

and the source root r that solves (j2^ is the quadrature rule of Nikolov m- 
The domain Cl C is generated as follows. For this quadrature rule we know 
the nodal pattern, e.g., for n odd, n = 2/c + 1, we have ( [2Q| ), therefore 

in 5 • • • 5 ) G [xo,Xi] X ••• X [Xk.Xk^i] X [Xk.Xk^i] X ••• X [Xn-l,Xn] (22) 

and for the weights we use (a very rough) range [0, b — a]. Combined together, 
the source domain is 

O = [xo, Xi] X • • • X [Xn-l,Xn] X [0, 6 — a] X • • • X [0, 6 — tt] . 

^-V-^ -V-^ (23) 

m m 


Remark 3. Eq. (23) gives a very loose bound on i = 1,..., m. However, this 
part of Cl is not affected by the nodal pattern and thus less important since it 
does not influence the change of the pattern of the polynomial system. A tighter 
bound could eventually serve as a better stopping criterion in cases, where the 
optimization is required to keep the root inside the domain (see Section 4.5) but 
we expect this would bring only marginal computational gains. 

There are several differences between the homotopy framework we propose 
and the classical setting mi: 


• There is only one root that we follow: the quadrature rule Q of More¬ 
over, we know the quadrature rules (the source rule, the intermediate ones 
and the target rule) are unique for geodesic knot transformations. This 
fact is a great advantage as there is no danger of numerical instabilities 
like jumping from one traced path to another as is the case in classical 
homotopy continuation methods. 

• The systems do not require linear blending as in [14]. In our case, the 
continuous transition between the systems is governed by the change of 
the knot vector. 
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Figure 6: Nodal pattern change. Top: evolution of the nodes (green) and knots (blue) is 
shown. At time instant ti, the source nodal pattern changes which corresponds to the fact 
that r 2 leaves {x 2 ^x^) by crossing the knot X 2 - Due to symmetry, T 3 leaves (^ 2 , 3 : 3 ) and the 
new nodal pattern is p = (1,1, 0,1,1). The new nodal pattern requires a change of the system 
( |21| ), e.g., T 2 starts affecting the basis function D 2 and the second equation of the system 
must be updated (right). Bottom: A zoom-in on the critical instant t = ti. 


• The domain O also changes in time because the range of every node is 
determined by the knot positions and these vary. 

• The system is only locally polynomial. When a node leaves its subinterval, 
the nodal pattern of the quadrature rule changes, a new polynomial system 
has to be built and solved. 

The last issue is important. The most difficult part of any optimal spline 
quadrature rule is to derive the correct layout of the nodes (the nodal pattern). 
Since the splines are piece-wise polynomials, the systems that we build and solve 
are also only “piece-wise polynomial”, i.e. locally polynomial and the main 
difficulty is to select the “right pieces”. However, homotopy continuation leads 
us to the correct nodal pattern automatically, by tracing down the trajectory 
of optimal quadrature rules as the knot pattern evolves. That is, the generated 
quadrature rules are exact (up to machine precision, see Section]^ and optimal 
for all the intermediate steps of the transformation. 

4 . 3 . Irregular nodal patterns 

We now discuss the changing the nodal patterns, see Fig.[^ This is the situ¬ 
ation when some (one or more) node(s) crosses the boundary of its subinterval. 
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That is, the node overlaps with a knot and the nodal pattern of the quadrature 
becomes irregular. We recall that is a hypercube in and the system 
(21) can be seen as an intersection problem of 2 m hypersurfaces inside 
which has at time t = 0 only one root r(0). During the continuation, as time 
evolves, 0(t) is changing and so does the root r(t) that is being traced. 

The situation when a node coincides with a knot corresponds to r(t) leaving 
D(t), i.e., there exist some ti G [0,1] such that r(ti) is on the boundary of the 
hypercube D(ti), r(ti) G D(ti). For simplicity, we assume that r(ti) lies on a 
hyperface of D(ti) which corresponds to the fact that only one node becomes a 
knot, ri{ti) = Xj{ti). Let the current sub-interval where Ti(t) lies in for t < be 
We know the hyper-face that is being crossed, Ti{ti) = Xj(ti), 
so the algorithm changes the nodal pattern by switching the sub-interval of 
from {xj-i^Xj) to {xj^Xjj^i). Then the system (21) is updated accordingly. 

If r(ti) lies on a hyper-edge, i.e., two (or more) nodes become knots at the 
same time, the situation is similar because for every node we know the current 
interval that the node leaves and the new one that it enters. 


4.4- Quadrature rule traeing 

The goal is to trace down a curve r(t) G t G [0, 1] knowing the initial 

point r(0), for example the quadrature rule of Nikolov [TT]. We know the source 
and the target knot sequences, ^ and (|^, respectively. As the number 
of non-degenerate intervals varies from n to A", we simplify the notation by 
omitting the subscripts and write T(0) := and T(l) = Tat. We recall that 
the number of interior knots and total knots remains constant. We further 
discretize time and build M — 1 intermediate knot vectors 

= = (24) 

and consequently discretize r(t) as a polyline {r*}^Q. We follow two different 
types of generation of the intermediate knot sequences, see Section The 
geodesic path, where each knot moves linearly in t to its target position, and 
along-the-edge paths, where only a single knot moves at each time instant, see 
Fig-i The first path is unique; the number of the paths of the second kind 
equals the number of “bottom-to-top” paths on a hypercube in (2n + 4)!. 

The latter paths generate spline spaces with mixed continuities (various knot 
multiplicities). 

Remark 4. The source quadrature rule is unique and so is the target rule. There¬ 
fore the target rule should get derived independently on the path between T(0) 
and T(l). However, to the best of our knowledge, there are no theoretical re¬ 
sults on the quadrature rules for splines with mixed continuities. The results of 
Micchelli and Pinkus [9] apply only for families of splines with uniform continu¬ 
ity. Our target quadrature rule gets derived independently of the path chosen, 
by parsing these spaces with mixed continuities as seen later in Section This 
result provides numerical evidence for the existence of optimal rules for these 
mixed continuity spaces. 
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Algorithm 1 


GaussianQuadrature([a, b],N) 


1: INPUT: compact interval [a, b] and odd number of uniform segments N 


2 

3 

4 

5 

6 

7 

8 


9 

10 

11 

12 

13 

14 

15 

16 

17 

18 

19 

20 


n=[f] + l; 
^0 


= the source rule Q on [a, b]; 
build = 0 from ([2^ over (23); 

build {T(ti)}^Q, Section ^ 
for i = 1 to M do 
build = 0 over 

solution of = 0 by multivariate Newton-Raphson 


opt • 

with initial guess r' 

if G then 


;-i. 


if 


< 5 then 

■■ opt "! 


^ opt 

r* := r 
else 

subdivide and 

M := M + 1; and go to line 7; 

end if 
else 

update and and go to line 7 

end if 
end for 

OUTPUT: r^, the Gaussian quadrature rule for S ^2 


/* the quadrature rule for X{ti) 
/* finer stepsize 

/* nodal pattern changes 


*/ 

*/ 


*/ 


4 . 5 . Implementation 

In our implementation, the time-discretization of the intermediate knot vec¬ 
tors is uniform. In the case of the geodesic knot transformation, we sample 
uniformly the diagonal (source-target) of the hypercube, see Fig.|^ In the case 
of the along-the-edge transformation, we sample uniformly every hyperedge. 

However, one can apply a non-uniform stepsize, e.g. this stepsize could get 
increased in cases where the nodes differ from the knots by more than a certain 
threshold. Reversely, a finer modification of the knot vector would be beneficial 
in cases when a node approaches a knot, i.e., in cases that precede a change of 
the nodal pattern. Such an analysis goes beyond the scope of this paper. We 
set M = 200 for the geodesic path and M = 20{N — 1) for along-the-edge paths. 

The tracing algorithm sequentially transforms the knot sequence PI{ti) and 
updates the system (21), = 0 considered over the domain Oh The exact 

quadrature rule from the previous iteration G is taken as an initial 
guess and a multivariate Newton-Raphson iterative scheme is applied to get the 
root If r^p^ lies in the current domain 0% the next knot sequence 
is taken with r* := r^p^ and the algorithm proceeds iteratively. When r^p^ ^ 0% 
this implies that the nodal pattern has changed as explained in Section [43] In 


these instances, the algorithm swaps to a different polynomial system defined 
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above . Then r* := and the Newton-Raphson optimization is applied 
again with the new system = 0 . If the optimization fails, i.e., if 

W^optW > ^ an additional intermediate knot sequence is inserted between 
and and the original point and the original system = 0 are 

considered again. Otherwise the nodal pattern is updated and := . The 

numerical threshold was set 5 = 10“^^. The whole procedure is summarized in 
Algorithm 


Remark 5. Algorithm requires the number of elements N to be odd. This 
is because the associated source space has all interior knots of multiplicity two 
and therefore N = n n — 1. This is in accordance with the result of Micchelli 
(17) that states the optimal rule for S ^2 guaranteed only for odd N. 


5. Numerical examples 


In this section, we show the results of our algorithm and derive Gaussian 
quadrature rules for various target spaces S ^2 • The evolution of the quadrature 
rules for = 11 is shown in Fig. Starting with the Gaussian quadrature rule 
for a spline space with uniformly distributed double knots as a source rule, the 
target rule for S ^2 traced by the homotopy continuation-based Algorithm 
deriving the optimal rules also for all the intermediate spaces. 

The evolution of the nodes is visualized as a set of planar curves, the evo¬ 
lution of weights is represented by the corresponding space trajectories. The 
error of the rule Q is measured in terms of the Euclidean norm of the vector of 


the residues of the system (21), normalized by the dimension of the system 


Ar+3 




(25) 


i=l 


Our results coincide with those of [2] (up to twelve digits precision shown 
therein), where the correct layout of nodes was conjectured as 


p = (l,1,0,1,0,1...,1,0,1,0,1,1), 


(26) 


that is, the first two boundary elements contain a single node each, whilst the 
middle elements follow a node-gap pattern. Our algorithm yielded the same 
nodal pattern, see Fig. Their algorithm is a trial-error scheme that chooses 
the first node and, under the assumption of the nodal pattern (26), computes 
iteratively the remaining nodes and weights. Based on the error of the rule, 
the first node is modified and the next iteration is invoked. Such an approach 
relies heavily on the initial node and, since the problem is highly non-linear. In 
our experience, there is no guarantee that the scheme of [2] will converge to a 
correct rule. We emphasize that our approach derives the layout of nodes (26) 


automatically and is general because it is not limited to only uniform target 
knot sequences. 
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Figure 7: Geodesic continuation of the Gaussian quadrature rules for N = 11. (a) The 

evolution of a Gaussian quadrature rule from the source space S'^ ^ with n = 6 is shown. As the 
knot vector changes (blue paths), the nodes (green) and the weights (red) move continuously 
from their source position (t = 0) to the target configuration (t = 1). The evolution of 
m = n + 1 = 7 Gaussian nodes and weights on [a, b] is shown from the top view (b) and 
the side view (c). (d) A zoom-in on the Gaussian quadrature rule for a specific time instant 
t = 0.25. The current nodes (green), weights (red) and knots (blue) are shown, (e) The target 
uniform knot configuration and the target rule at t = 1. 
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Figure 8: The evolution of Gaussian quadrature rules ( |19t for various source N are shown. 
Top: The evolution of the nodes (green), weights (red) for geodesic knot transformation (blue). 
Middle: The evolution of nodes in time and the target rules for the time instant t = 1. Bottom: 
The rule for N = 39; observe the convergence to the midpoint rule of Hughes et al. [8], cf. 
Table [T] 


Also observe the convergence to the midpoint rule of Hughes et al. [8], see 
Table where the limit weight for N = 39 {n = 20) is 

^ = 0.051282 (27) 

Oc/ 

The data in Tableare shown with double-precision (16 decimal digits). Ob¬ 
serve the slow convergence, e.g., for N = 39, the results meet the half-point-rule 
limit with only five digits of accuracy. This fact is not a limitation of our work, 
on the contrary, it shows that the midpoint rule is an approximation of our 
exact rule and can be used only when N is large enough. 

The evolution of the Gaussian nodes and weights for various N is shown in 
Fig-i The modification of the knots is geodesic, see Section i.e., every knot 
transforms to its target destination along the shortest possible path. As previ¬ 
ously discussed in Section the source and the target rules are both unique. 
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Table 1: Nodes and weights for Gaussian quadrature rules ( |19| > with double-precision 
for uniform knot distribution for various target N are shown. To compare the results 
with algorithm of Nikolov p], the interval was set as [a,b] = [0,1]. Due to the 
symmetry, only the first [ ] + 1 nodes and weights are displayed. The error 

||r|| of the rule is measured as the Euclidean no rm of the vector of the residues, 
normalized by the system’s dimension A/” + 3, see ( |25| >. 


i 

Ti 

OJi 

l|r|| 


N : 

= 3 


1 

0.1086264370680297 

0.2720231005023455 

7.90-2“ 

2 

0.5 

0.4559537989953090 



N : 

= 5 


1 

0.0669578918742195 

0.1698605936669416 

1.04-1® 

2 

0.3275898516368645 

0.3301394063330584 



N : 

= 7 


1 

0.0479188107803577 

0.1216810800700958 


2 

0.2358921494969001 

0.2408185184939348 

1 . 95 - 1 ® 

3 

0.5 

0.2750008028719389 



N : 

= 9 


1 

0.0372757529111283 

0.0946622477445919 


2 

0.1835904624135774 

0.1876252194189693 

2.08-1® 

3 

0.3904233866079767 

0.2177125328364388 



N = 

= 11 


1 

0.0304987043023585 

0.0774523185174377 


2 

0.1502181009517147 

0.1535325192913209 

6 . 68 - 1 ® 

3 

0.3195393932155687 

0.1783894870783702 


4 

0.5 

0.1812513502257421 





= 39 


1 

0.0086022074347388 

0.0218455595269063 


2 

0.0423693959303822 

0.0433045545577068 


3 

0.0901289847662636 

0.0503213631747089 


4 

0.1410569521267253 

0.0512021143533085 


5 

0.1923101843694322 

0.0512756766459810 


6 

0.2435899416018961 

0.0512815446928528 

1.02-^’’ 

7 

0.2948718106031808 

0.0512820110347811 


8 

0.3461538474036372 

0.0512820480845737 


9 

0.3974358975351839 

0.0512820510280155 


10 

0.4487179487257872 

0.0512820512617426 


11 

0.5 

0.0512820512788446 
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t = 0 





(4,1,3, 2) 





1 


Figure 9: Path in depe ndence. The evolution of the Gaussian 
quadrature rules ( |19| > for various modifications of the knot 
vector for N = 5 is shown. At every time instant, only one 
knot changes its position. Five different evolutions of the 
Gaussian nodes are shown (green). Starting with two double 
knots, the quadruplet encodes the particular order of the four 
knots to be moved towards their target position (uniform 
spacing at t = 1). The 3D figure compares the (4,1,3,2)- 
evolution of the weights (blue) with the geodesic evolution 
(yellow), cf. Fig.[^ All paths derive the same target rule {t = 
1) visualized by red dots (weights) and green dots (nodes). 


Therefore, any path deforming the uniform double knot sequence Q into 
the uniform knot sequence Tat (§ must produce the same rule. The results of 
five random along-the-edge knot modifications are shown in Fig. Contrary 
to the geodesic knot deformation, only one knot moves while the other stay 
fixed during the evolution. Depending on the order of how the knots are moved, 
there exist {N + 5)! possible paths. The results of the target rule for various 
along-the-edge knot vector deformations are shown in Table The results are 
independent on the path which validates the correctness of our algorithm. 

An example of Gaussian rules for non-uniform target knot sequences is shown 
in Fig. 10 The space is kept (7^, i.e., all the knots keep their multiplicity two. 


The target knot sequence was required to have larger elements close to the 
boundary and smaller ones in the interval’s center. Unlike our previous result 
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Table 2: Path independence. The target quadrature rule obtained by different along- 
the-edge paths of the knot vector. The results for five paths shown in Fig. are 
displayed with twenty digits accuracy. The target nodes and weights match up to 
eighteen digits. 



Tl, T2 

LUi, LU2 

(1,2,3,4) 

0.06695789187421950918 

0.32758985163686446374 

0.16986059366694164265 

0.33013940633305835725 

(1,2,4,3) 

0.06695789187421950918 

0.32758985163686446374 

0.16986059366694164265 

0.33013940633305835725 

(1,4, 2,3) 

0.06695789187421950918 

0.32758985163686446374 

0.16986059366694164265 

0.33013940633305835725 

(2, 3,4,1) 

0.06695789187421950915 

0.32758985163686446341 

0.16986059366694164265 

0.33013940633305835787 

(4,1,3, 2) 

0.06695789187421950915 

0.32758985163686446342 

0.16986059366694164265 

0.33013940633305835721 


[T] where finer spacing close to the boundary guaranteed the same nodal pattern 
for the whole family of symmetrically stretched knot sequences, here we can 
observe a change of the nodal pattern of the rule. 

Throughout the paper, we used the rule of Nikolov m as our source rule. 
However, one may start, e.g., with the classical Gaussian rule for polynomials, 
see Fig. In such a case, the knots have multiplicity four and one needs to 
count the proper number of initial elements to get the desired number of target 
elements N. This counting obeys the rule of Micchelli [ 9 ], see ( p! 7 [ ). The results 
show the stability of our Algorithmic The results from Fig.[^both correspond 
to values shown in Table [C for N = 9. Both target rules derived from different 
source rules match up to sixteen decimal digits. 


6. Conclusion 


We have introduced a new methodology to compute Gaussian quadrature 
rules. Starting with a known Gaussian quadrature rule for a specific spline 
space, the rule is interpreted as a point in a higher dimensional space, a zero 
of a specific polynomial system. A homotopy continuation-based algorithm has 
been presented that numerically traces the Gaussian quadrature rule as the 
knot vector, and consequently the whole spline space, is continuously modified. 
We have recovered the Gaussian quadrature rule of Nikolov [2] for the uniform 
cubic spline spaces, a rule that was derived under an assumption of the 


conjectured nodal pattern (26). 


The continuation approach shows the connection between different Gaussian 
rules. This concept eliminates the need for the construction of particular rules 
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Figure 10: The evolution of the Gaussian rules for spaces with non-uniform target knot 
sequences. Five interior source double knots (blue) are moved towards the interval’s center. 
The parameter q is the shrinking ratio; it is the ratio of lengths of two neighboring elements 
when moving towards the center. Top: The evolution of the rule in time visualized in 3D: 
the trajectories of the weights (red), nodes (green) and knots (blue) are shown. Middle: The 
evolution of the nodes. For all q shown, T 2 starts (t = 0) in the second elements but ends 
(t = 1) in the first element. Bottom: The target Gaussian rules for particular non-uniform 
knot sequences. 

for different kinds of non-uniform knot vectors. Thus, work such as mm is 
hereafter needed only to validate numerical results obtained with this method. 
We plan to further exploit the proposed technique to derive the Gaussian rules 
for the particular spline spaces of a high interest such as 84^0 with non-uniform 
knot sequences, or S'e,!, spaces frequently appearing in the Galerkin method 
when building the mass and stiffness matrices. 
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